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We present a novel quantum-classical approach to non-adiabatic dynamics, deduced from the cou¬ 
pled electronic and nuclear equations in the framework of the exact factorization of the electron- 
nuclear wave function. The method is based on the quasi-classical interpretation of the nuclear wave 
function, whose phase is related to the classical momentum and whose density is represented in terms 
of classical trajectories. In this approximation, electronic decoherence is naturally induced as effect 
of the coupling to the nuclei and correctly reproduces the expected quantum behaviour. Moreover, 
the splitting of the nuclear wave packet is captured as consequence of the correct approximation 
of the time-dependent potential of the theory. This new approach offers a clear improvement over 
Ehrenfest-like dynamics. The theoretical derivation presented in the Letter is supported by numerical 
results that are compared to quantum mechanical calculations. 
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The theoretical description of phenomena such as vi¬ 
sion [1], photo-synthesis [2], photo-voltaic processes [3], 
proton-transfer and hydrogen sforage [4] is among 
fhe mosf challenging problems in Condensed Matter 
Physics and Theoretical Chemisfry. The underlying 
quanfum dynamics of elecfrons and nuclei exhibif a 
non-adiabafic characfer, meaning fhaf if cannof be ex¬ 
plained by employing fhe Bom-Oppenheimer (BO) ap- 
proximafion. In fhis respecf, fhe major challenge for fhe- 
ory resides in fhe explicif freafmenf of elecfronic excifed- 
sfafe dymamics coupled fo fhe nuclear motion. While 
mefhods fhaf refain quanfum feafures of fhe nuclear dy¬ 
namics [5] are fhe mosf accurafe fo address fhis prob¬ 
lem, fhey cannof be applied fo sysfems wifh hundreds, 
or even fhousands, of afoms. Therefore, a freafmenf 
of nuclear dymamics in ferms of (semi)classical frajec- 
fories [5-9] represenf fhe mosf promising and numeri¬ 
cally feasible approach for acfual calculafions. Despife 
fhe greaf efforf fhaf has been devofed over fhe years fo 
fhe developmenf of such mefhods, acfual applications 
are still limifed [10]. Well-known issues are connecfed fo 
fhe lack of, or incorrecf accounf for, decoherence and fo 
fhe inabilify of reproducing fhe spatial splitting of a nu¬ 
clear wave packef, as in Ehrenfesf-like dymamics. In fhe 
sfudy of elecfronic non-adiabafic processes, fhese prob¬ 
lems can resulf in wrong predicfions for quanfum pop¬ 
ulations and in unphysical oufcomes for fhe nuclear dy¬ 
namics. 

We have recenfly proposed a new formalism fhaf can 
be employed fo overcome fhe above issues, fhe so-called 
exacf factorization of fhe elecfron-nuclear wave func- 
fion [11]. In fhis framework, fhe full wave function 
is written as fhe producf of a nuclear wave function 
and an elecfronic factor wifh a paramefric dependence 
on fhe nuclear configuration. Coupled equations drive 


fhe dymamics of fhe fwo componenfs of fhe wave func- 
fion. In particular, a fime-dependenf Schrodinger equa¬ 
tion (TDSE) describes fhe evolufion of fhe nuclear wave 
function where fhe effecf of fhe elecfrons, beyond BO, 
is accounted for in a single, fime-dependenf, pofenfial. 
Compared fo a formulafion in ferms of multiple sfafic 
adiabatic (or BO) pofenfial energy surfaces (PESs), fhe 
advanfage of fhis formulafion is evidenf: when fhe clas¬ 
sical approximation is infroduced, fhe force driving fhe 
nuclear evolufion can be uniquely determined from fhe 
gradienf of fhis fime-dependenf pofenfial. [12] 

In previous work we have: (i) analyzed fhe feafures 
of fhe fime-dependenf pofenfial [13] in fhe confexf of 
non-adiabafic profon-coupled-elecfron-fransfer, in or¬ 
der fo pinpoinf fhe properties fhaf need fo be accounfed 
for when infroducing approximafions; (ii) defermined 
fhe suifabilify of fhe classical and quasi-classical freaf¬ 
menf [14] of nuclear dynamics, in a sifuafion where fhe 
elecfronic effecf can be faken info accounf exacfly; (iii) 
derived an independenf-frajecfory (IT) mixed quanfum- 
classical (MQC) algorifhm [7, 8] fo solve fhe coupled 
elecfronic and nuclear equations (from fhe factorization) 
in a fully approximafe way. In particular, fhe IT-MQC 
scheme has been obfained as fhe lowest-order approx¬ 
imation, in an expansion in powers of h of fhe nu¬ 
clear wave function in fhe complex-phase represenfa- 
fion. Eurfher investigation [15] has shown, however, 
fhaf correcfions are required if fhe nuclei exhibif a quan¬ 
fum behavior relafed fo a non-adiabafic evenf, e.g. fhe 
spliffing of a nuclear wave packef after fhe passage 
fhrough an avoided crossing. 

The aim of fhis Letter is fo go beyond fhe IT-MQC 
algorifhm of Ref. [7, 8]. We have derived a coupled- 
frajecfory (CT) MQC algorifhm able fo reproduce fhe 
feafures of fhe fime-dependenf pofenfial, by evolving an 



ensemble of classical trajectories to mimic the quantum 
evolution of the nuclei. Electronic populations, deco¬ 
herence and spatial splitting of the nuclear wave packet 
are correctly reproduced when the new scheme is em¬ 
ployed, as will be demonstrated below. 

The exact factorization approach consists in writing 
the solution, 'k(r, R, t), of the TDSE = ztidt'k, as 
the single product 'I'(r, R, t) = $R(r, t)x(R, t), where 
is 31^ electronic factor parametrically depend¬ 
ing on the nuclear positions and x(R, t) is a nuclear 
wave function. Here, H = Tn + Hbo the Hamilto¬ 
nian describing the system of interacting electrons and 
nuclei, with the nuclear kinetic energy and Hbo the 
BO Hamiltonian containing all interactions among the 
particles and the electronic kinetic energy. The posi¬ 
tions of Ne electrons and iV„ nuclei are represented by 
the symbols r and R, respectively. The product-form of 
T* is unique, up to within an (R, t)-dependent gauge¬ 
like transformation, if the partial normalization condi¬ 
tion, f dr|<I)R(r, = 1 VR, t, is imposed. The evolu¬ 
tion of the two components of the full wave function is 
governed by an electronic equation. 


HbO + l^en[^R, X] 


e(R,t) 


<t>R — ihdt^R, 


( 1 ) 


and a nuclear equation. 
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X = ihdtx, (2) 


which are exactly equivalent to the full TDSE. The 
electron-nuclear coupling operator. 
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of R^^^(t),t. Nuclear d 5 mamics will be sampled us¬ 
ing trajectories, meaning that we track the evolution 
of a nuclear wave packet by looking at how the tra¬ 
jectories evolve in time. Information about the nuclear 
space R is only available at the instantaneous positions 
along the classical paths. It follows that we will not 
be able to calculate partial time derivatives, but only 
total time derivatives, by using the chain rule d/dt = 
with = r[,^^ (t) the nuclear veloc¬ 
ity. Henceforth, the superscript (/) will be used to indi¬ 
cate a spatial dependence, e.g., = A„(R‘''^^(t), t). 

The main steps in the derivation of the new CT-MQC 
scheme are the following: (a) we approximate the TD- 
PES, to avoid expensive calculations of second-order 
derivatives of the electronic wave function with-respect- 
to the nuclear coordinates; (b) we fix the gauge free¬ 
dom; (c) we introduce a quasi-classical interpretation of 
the nuclear wave function, whose phase is connected to 
the classical momentum and modulus reconstructed in 
terms of Gaussian wave packets; (d) we expand the elec¬ 
tronic wave function on the adiabatic basis (Bom-Huang 
expansion), exp[(z/fi) 7 f^^(t)]p[^^ 

hence a set of partial differential equations for the coef¬ 
ficients of the expansion will be coupled to the nuclear 
equation. Also the full wave function can be expanded 
on the adiabatic basis, with coefficients T/(R, t), for the 

exact expression, or for the quantum-classical 

case, and will be referred to as BO-projected wave pack¬ 
ets. 

(a) In the expression of the TDPES we neglect the 
contribution of ($R(t)|17en|‘l’R(t))r- Notice that the ex¬ 
pectation value on <I>r of the second line of Eq. (3) is 
zero by construction, thus the neglected term in the ex¬ 
pression of the TDPES contains the second-order varia¬ 
tions of the electronic state with-respect-to the nuclear 
coordinates, which is small [16-18] compared to the 
first-order. Therefore, the TDPES is approximated as 
e(R, t) ~ eo(R, t) + eTD(R, t) and Eqs. (1) and (2) be¬ 
come 


represents the effect of the nuclei on electronic dynam¬ 
ics; in turn, the time-dependent vector potential, 

A„(R,t) = ($g(t)| - 

and time-dependent PES (TDPES), 


account for the electronic back-reaction on the nuclei 
in a Schrodinger-like equation. These potentials are 
uniquely determined [11] up to within a gauge trans¬ 
formation. 

The CT-MQC scheme adopts a description of nuclear 
d 5 mamics in terms of classical trajectories, R*^^^ it), thus 
all quantities depending on R, t will become functions 


ih^d) = ^ ^(6) 
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respectively, where the symbol is used to indicate 
the full time-derivative of the electronic wave function 
and will be specified below. The equations have 
been cast in such a way that the first terms on the 
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right-hand-side are exactly the same as in the Ehren- 
fest scheme [19]. The additional terms are corrections, 
whose effect will be now investigated. 

(b) In deriving these expressions for fhe evolution of 

fhe elecfronic wave function, Eq. (6), and for fhe classical 
nuclear force, Eq. (7), fhe gauge freedom has fbeen fixed 
by imposing -I- €^^(t) + = 0- 

(c) The corrections beyond-Ehrenfesf in 

Eqs. (6) and (7) confain a ferm — 

which we will refer fo as 
quantum momentum. The reason for fhis choice lies in 
fhe following expression 


-ihv^xiRyt) 

x(E,i) 


A„(E, t) 


p(J)(t) + r(J\t) 


( 8 ) 


for fhe ferm in Eq. (3) fhaf explicifly depends on fhe nu¬ 
clear wave function. Such ferm has fo be approximafed 
when a frajecfory-based freafmenf is adopfed. In facf, 
Eq. (8) has been obfained by writing fhe nuclear wave 
function in polar form, x = ^rid fhen idenfi- 

fying (quasi-classically) -\- A„ = P„, wifh fhe 
classical nuclear momenfum and P„ = F„ from Eq. (7). 

(d) If compared fo fhe Ehrenfesf scheme, fhe imple- 
menfafion of fhe CT-MQC algorifhm based on Eqs. (6) 
and (7) requires only two additional steps: the cal¬ 
culation of (i) and (ii) We employ fhe 

Bom-Huang expansion of in order fo express 

fhe ferm (i) using fhe derivafives of fhe expansion co- 
efficienfs, indicafed by fhe symbols c\^\t), and fhe 
non-adiabafic coupling vecfors. The approximafion 
~ {i/h)W used here, wifh Xi^\t) fhe 
phase of c[^\t), is consisfenf wifh previous analy¬ 
sis reporfed in Refs. [13, 14]. Moreover, based on 
quasi-classical considerations described in defail in 
fhe Supplemenfal Maferial, we furfher approximafe 
~ drV The ferm (ii) is calculafed 
assuming fhaf fhe nuclear densify is a combinafion of 
Gaussian-shaped wave packefs, each corresponding fo 
a given adiabafic sfafe. Notice fhaf fhis approximafion 
is nof used in general in fhe algorifhm, buf only fo esfi- 
mafe fhe quanfum momenfum. Eor a fwo-sfafe model, 
becomes [15] a linear funcfion in fhe region where 
= \cl^\t)\'^ ^ 0,1, while if is sef fo zero else¬ 
where (see fhe Supplemenfal Maferial and fhe discus¬ 
sion below). The generalization of fhis approximafion fo 
multiple sfafes is sfraighfforward and will be presenfed 
elsewhere [20]. The paramefers of such linear funcfion 
are fhe slope and fhe j/-infercepf, where fhe former is de- 
fermined analytically by using Gaussian-shaped nuclear 
wave packefs and fhe laffer is obfained by enforcing (fhe 
reasonably physical condition) fhaf no population ex¬ 
change occurs when fhe non-adiabafic coupling vecfors 
are zero. Information abouf fhe positions of all frajec- 
fories af a given fime is required when evaluating fhese 


fwo paramefers, fhus resulting in a procedure beyond 
fhe IT-MQG approach: fhe classical frajecfories carmof 
be evolved independenfly from each ofher, fhey are cou¬ 
pled. 

The major advanfage of fhe GT-MQG scheme devel¬ 
oped here is fhaf fhis procedure nafurally incorporafes 
decoherence effecfs. In fhe following we shall discuss fhis 
feafure in defail. Affer fhe nuclear wave packef has leff a 
region of sfrong non-adiabafic coupling, fhe population 
of ff*o /-fh BO sfafe changes in fime 
as 



In fhis region, fhe expression of fhe vecfor pofenfial re¬ 
duces fo A^v \t) = J2i since fhe non- 

adiabafic coupling vecfors are negligible. In Eq. (9) we 
observe fhaf, once (f) has approached fhe values 0 
or 1, fhe ferm on fhe righf-hand-side becomes zero, fhus 
fhe elecfronic population remains consfanf (fo 0 or 1) 
V 1. This is a clear indication of decoherence, since fhe 
(squared-modulus of fhe) off-diagonal elemenfs of fhe 
elecfronic densify mafrix, offen used as a measure of 
elecfronic coherence, become zero. Therefore, fhe cor- 
recfion ferms beyond-Ehrenfesf in Eqs. (6) and (7), pro- 
porfional fo fhe quanfum momenfum, will be referred fo 
as decoherence terms. Obfaining fhis feafure is a clear 
improvemenf over fhe Ehrenfesf approach and, like¬ 
wise, over fhe IT-MQG approach [7, 8] deduced from 
fhe exacf facforizafion. Decoherence nafurally appears 
by including dominanf correcfions in fhe expression of 
fhe nuclear wave funcfion, leading fo fhe appearance of 
fhe quanfum momenfum. 

Numerical resulfs obfained by implementing fhe 
above-described mefhod are shown below in compari¬ 
son fo exacf calculafions. We discuss fhe performance 
of fhe GT-MQG algorifhm in comparison fo Ehrenfesf 
d 5 mamics for a fwo-sfafe problem [21] involving fhe 
passage of fhe nuclear wave packef fhrough a single 
avoided crossing, case (1), and in comparison fo frajec- 
fory surface hopping (TSH) for a fwo-sfafe problem in¬ 
volving fhe reflecfion of fhe nuclear wave packef from 
a pofenfial barrier and ifs consequenf spafial spliffing, 
case (2), commonly known as Tully-3 model [9]. The 
compufafional defails for fhese fwo problems are given 
in fhe Supplemenfal Maferial. 

The TDPES for model case (1) is shown in Eig. 1 (up¬ 
per panels). If develops sfeps and fhe nuclear wave 
packef correcfly splifs af fhe avoided crossing (Eig. 1, 
lower panels). If is worfh nofing fhaf Ehrenfesf d 5 mam- 
ics complefely misses fhe spliffing, as we have shown 
in Ref. [8]. Eurfhermore, despife fhe facf fhaf Ehrenfesf 
d 5 mamics properly reproduces fhe populafions of fhe 
elecfronic sfafes, as shown in Eig. 3 (upper leff panel), 
if does nof capfure decoherence. Qn fhe confrary fhe 
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CT-MQC procedure slightly underestimates the non- 
adiabatic population exchange but correctly reproduce 
decoherence (Fig. 3, lower left panel). 

In Fig. 3, we have used the quantity 

measure of de¬ 
coherence, whose quanfum equivalenf is 
/ dR;Pi(R, f)p 2 (R, f)|x(S, Here, fhe nuclear 

densify has been replaced by ifs "classicaF' approxima- 
fion,i.e. |x(a,f)P ~ iV-h 5(a - 

We show in Fig. 2 (lower panels) fhaf fhe CT-MQC 
algorifhm reproduces fhe spliffing of fhe nuclear wave 
packef due fo fhe reflecfion from fhe barrier. In facf, fhe 
TDPES develops fhe well-sfudied [13] sfeps (Fig. 2, up¬ 
per panels) fhaf bridge piecewise adiabatic shapes and 
allow frajecfories in differenf regions of space fo feel 
differenf forces. This feafure is fhe sfrengfh of a pro¬ 
cedure based on fhe exacf facforizafion: a single time- 
dependenf pofenfial generafing very differenf forces in 
differenf regions of space. If is known [9] fhaf TSH as 
well is able fo capfure fhe reflecfion evenf for a low ini¬ 
tial momenfum of fhe nuclear wave packef, buf suffers 
from over-coherence [10, 22]. In facf, as shown in Fig. 3, 
TSH complefely misses decoherence, whereas fhe CT- 
MQC scheme nof only reproduces fhe populations of 
fhe elecfronic sfafes as functions of fime (Fig. 3, upper 
righf panel), buf can also capfure elecfronic decoherence 
(Fig. 3, lower righf panel). The comparison between 
exacf and CT-MQC resulfs is overall remarkable and a 
clear sfep forward in comparison fo ofher mefhods. 

In fhis Leffer we have proposed a CT-MQC scheme 
based on fhe exacf facforizafion formalism and fesfed if 
on a f 5 q)ical example of elecfronic non-adiabafic process. 
The resulting equations give additional ferms compared 
fo Ehrefensf d 5 mamics, fhaf appear fo be responsible for 
decoherence. The comparison of fhe CT-MQC scheme 
wifh full quanfum mechanical resulfs shows fhaf we can 
correcfly predicf bofh elecfronic and nuclear properties: 
population d 5 mamics, nuclear wave packef spliffing and 
decoherence. Non-adiabafic fransifions are induced by 
the classical nuclear momentum, the zero-th order term 
of fhe fi-expansion of fhe nuclear wave function, and de¬ 
coherence is fhe effecf of fhe dominanf correcfions fo fhe 
momenfum. In addition, we have proven fhaf, as dis¬ 
cussed in our previous work [13,14], being able fo cafch 
fhe main feafures of fhe fime-dependenf pofenfial in an 
approximafe scheme resulfs in fhe correcf descripfion of 
fhe nuclear d 5 mamics. The major advanfages of our CT- 
MQC algorifhm over commonly used mefhods are: (1) 
fhe working equafions are concepfually and compufa- 
fionally as simple as Ehrenfesf equafions, and (2) a small 
number of frajecfories is required, because only inifial 
conditions are fo be sampled (no sfochasfic elemenf is 
introduced). Working in the framework of fhe exacf fac¬ 
forizafion allows fo sysfemafically improve previous ap- 
proximafions, as we have shown in fhis Leffer in com- 
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FIG. 1. Snapshots at different times for model case (1) of: 
(top) the gauge-invariant (GI) part of the scalar potential 
CG/(S, t) = {^R{t)\HBO + Uen\^R(t))r (Imes) Compared to 
its approximation (dots) used in the CT-MQC calcula¬ 

tions, in Hartree eh (the BO surface are plotted for reference as 
dashed black lines); (bottom) nuclear density |x(B:> t)P (black 
line) and BO-projected densities (dashed lines) compared to 
the values of the BO-projected densities along the trajectories 
(dots). 


t = 47.4 fs t = 68.9 fs t = 98.2 fs 



R (ao) R (ao) R (ao) 


FIG. 2. Same as Fig. 1 for model case (2). 

parison fo fhe IT-MQC of Ref. [7,8]. Along similar lines, 
fufure work will focus on including quanfum nuclear ef- 
fecfs, such as inferference, adopting a semiclassical rep- 
resenfafion of nuclear d 5 mamics. 

Partial supporf from fhe Deufsche Eorschungsge- 
meinschaff (SEE 762) and from fhe European Commis- 
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time (fs) time (fs) 


FIG. 3. Populations of the BO states (upper panels) and in¬ 
dicator of decoherence (lower panels) as functions of time for 
model case (1) (left) and model case (2) (right). 
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In the Letter we state that the main evolution equations of the coupled-trajectory 
mixed quantum-classical (CT-MQC) scheme can be cast such that corrections terms to 
Ehrenfest-like equations arise naturally when the quantum momentum is included. 
Here, we will present this result by employing the Born-Huang expansion to represent 
the electronic wave function <l>R(r, t) as a superposition of adiabatic states. Alternative 
expressions of Eqs. (6) and (7) of the Letter will be derived. Notice that this operation is 
simply a way to re-write Eqs. (6) and (7), no additional approximations are considered, 
rather than those already discussed in the Letter. 


EQUATIONS EMPLOYING THE BORN-HUANG EXPANSION 


The adiabatic states are the eigenvectors of the Bom-Oppenheimer (BO) Hamiltonian, 
Hbo> and are indicted by the symbols {v^r where Ngt is the number of states 

that will be included in the expansion. The corresponding eigenvalues are labeled as 
e^Q(R). The Born-Huang expansion is 




(l, t) = X] IC'i (R, 0 I exp 


i=i 


jpi (R,f) 


(i) 


(SM.l) 


where the coefficients G(R, t) have been written in terms of their moduli, |Oz(R, t) |, and 
phases, 7 /(R, t). 

The evolution equations for the coefficients C'z(R, f) and the classical nuclear force, 
from Eqs. (6) and (7), can be written now as 


Nst 

c\y) ■ Ecf «)d 


{!) 

v.lk 




and 






Nst 




(fc),(/) _ M) 

^BO ^BO j ^vlk 


^ (7) 

\x(n(t)\ 


(SM.2) 

(SM.3) 


(SM.4) 

(SM.5) 


with = {i/h)'P\j\t). Here, in order to adopt the same notation used in 

the main text, we indicate the R-dependence by introducing the label (/), the trajectory, 
and dropping the label R. New symbols have been introduced, namely 




and = 


(SM.6) 
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the standard non-adiabatic coupling vectors, defined as 










(SM.7) 


the gradient of the phases, of the expansion coefficients in Eq. (SM.l) 
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and the matrix 
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A (A»'(i) _ , (Sm.9) 

l=l 

In this last expression, the symbol A stands for a tensor product. Moreover, notice that, 
using the Born-Huang expansion in Eq. (SM.l), also the vector potential can be (exactly) 
re-written as 
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Nst 
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(SM.IO) 


where real electronic BO states have been assumed. 

Eqs. (SM.2) and (SM.4) are exactly the expressions appearing in the Ehrenfest algo¬ 
rithm, when the Born-Huang expansion is used to represent the electronic wave function. 
The additional terms depend on the quantum momentum, defined as 

n'Ht) = (SM.li) 


in Eq. (8) of the Letter. 

Eqs. (SM.2) to (SM.5) are the basic equations of the CT-MQC algorithm derived and 
discussed in the Letter. 


APPROXIMATION FOR THE GRADIENT OF THE EXPANSION COEFFICIENTS 

The quantity (t) appears in the evolution equations as consequence of the following 

approximations. The gradient of the expansion coefficients in Eq. (SM.l), namely 


V„c<''(i) = 






cP(t) 


cP(t) 


(SM.12) 


is approximated as 


V„C,<'>(() ~ jCP(t)V.T!‘’(t), 


p). 


(SM.13) 
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since we have already observed, based on full quantum calculations [1-3], that the contri¬ 
bution Vl 


Ci^{t) I is considerably different from zero only in regions where the nuclear 
density is small and thus where only few trajectories are expected. The error introduced 
this way can be systematically reduced by increasing the number of trajectories. 


In the calculations, we use an approximation to the term (t) based on the follow¬ 

ing quasi-classical considerations. When the Born-Huang expansion is also employed for 
the full wave function T (r, R, t), namely 


Nst 




R, 11 exp 


1=1 


J.SI {R,t) 


(i) 


(SM.14) 


the phases of the coefficients in Eq. (SM.l) and of the coefficients Fi (R, t) are related via 
the expression 


Fit) = 4'\t) - (SM.15) 

Here, once again, we have replaced the spatial dependence by the dependence on the 
trajectory I and we have used the symbol to indicate the phase of the nuclear 

wave function, as done in the Letter. Eq. (SM.15) straightforwardly follows from the 
factorization. It we identify the nuclear momentum calculated along the J-th trajectory 
as + Au\t) = Pu\t) and similarly as the momentum 

associated to the propagation of the /-th "BO-projected" wave packet (t) p on the /-th 
BO adiabatic surface, we can write 

Vy'V) ^ ?<"'''>(«) - (Pi,'>(i) - A<'>(i)) . (SM.16) 

Notice that is identified [1-3] as a "BO-projected" wave packet since from the 

factorization it follows that |xp = l-Hp. Taking the time-derivative of both sides of 
Eq. (SM.16) 


(SM.17) 

since on the right-hand-side we identify the forces associated to the motions of the "BO- 
projected" wave packet and of the nuclear wave function. The time-derivative of the term 
in parenthesis in Eq. (SM.16) is identically zero because in the chosen gauge the scalar 
potential is zero, thus all effects of the coupling to the electrons in the nuclear equation is 
contained in the vector potential. Indeed, Eq. (SM.17) is valid only in the gauge we chose 
to derive the CT-MQC algorithm. The expression for {t) used in our algorithm thus 

follows, 

F'(t) = VFF) = - /’ ‘IrVj'F- (SM.18) 
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APPROXIMATION FOR THE QUANTUM MOMENTUM 


We recall here the expression of the quantum momentum 




\xm\ 


in -!-^ 


(SM.19) 


which is written in the second equality in terms of the nuclear density In or¬ 
der to derive a practical way to estimate this expression, we assume that can 

be written as a sum of Gaussian-shaped wave packefs. In the simple case of fwo adia¬ 
batic states, only two Gaussian functions are considered, as schematically illustrated in 
Fig. SM.l. We underline again (as in the Letter) that this hypothesis is introduced only 
to estimate the quantum momentum and not to actually solve the classical nuclear equa¬ 
tion. We refer fo [4], where we show thaf the approximation we use here for the quantum 
momentum is consistent also with a representation of the nuclear wave packet in terms 
of coherent states. Two scenarios may occur: (1) the BO-projected wave packets overlap 



Supplemental Material, Figure SM.l. Schematic plot of the shape of the quantum momentum for 
various cases. The blue and red Gaussians represent BO-projected nuclear densities while the 
blue dashed lines are the quantum momentum (apart from the —ih multiplicative factor) corre¬ 
sponding to each configuration of the Gaussians. The grey areas highlight the regions where the 
nuclear density is considerably different from zero. 


and classical frajecfories may be located in the overlap region; (2) the BO-projected wave 
packets are well-separated, thus the overlap region is poorly sampled by the trajectories 
since the nuclear density is small. In case (1), we use a linear function to approximate the 
quantum momentum, as shown in the panels (A) and (B) of Fig. SM.l, within the over¬ 
lap region, while we set it equal to zero outside this region. Since outside the overlap 
region there is a low probability of finding frajecfories, as the nuclear density is small, 
the final resulf is nof strongly affected by this approximation. Moreover, as discussed 
in Eq. (9) in the Letter, the term in the electronic equation containing the quantum mo¬ 
mentum vanishes outside the overlap region, if in this region the non-adiabatic coupling 
vectors are zero. The inverse situation occurs in case (2), where the linear approximation 
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does not strongly affects the trajectories, because the overlap region is poorly sampled. 
This is shown in the panel (C) of Fig. SM.l. Outside this region, where the majority of 
trajectories is located, the quantum momentum is set to zero. However, this situation in 
general takes place after the nuclear packet has left the region of non-adiabatic coupling, 
where Eq. (SM.3), containing the quantum momentum, is identically zero as discussed in 
the Letter in Eq. (9) (this feature is related to the fact that if \Cj'^\t)\ = 1 all other 
\/k ^ I are identically zero). 


IMPLEMENTING THE CT-MQC ALGORITHM 

We describe the steps necessary to implement the CT-MQC: 

1. get the BO potentials along with and the non-adiabatic coupling 

vectors 

2. compute using the BO forces; 

3. collect information about the positions of all trajectories, (t) (this step implies 
that the trajectories caimot be propagated independently); 

4. calculate the quantum momentum, = a [Rd)(t) - 

where we use the symbol Ruit) for the crossing point between the two Gaussians 
(the y-tntercept), as illustrated in the Letter, namely by enforcing (the reasonably 
physical condition) that no population exchange occurs when the non-adiabatic 
coupling vectors are zero; 

5. evolve the coefficients {t) and the trajectories according to Eqs. (SM.2) to (SM.5); 

6. go back to step 1. 

COMPUTATIONAL DETAILS 

Test case (1) 

The model system for non-adiabatic charge transfer [5] analyzed in previous work [1- 
4, 6] has been also employed in the present study, as it exhibits the typical features of 
non-adiabatic processes and is exactly solvable, by integrating the full TDSE using the 
split-operator technique [7]. It thus provides a benchmark for the CT-MQC scheme. 
The system consists of three ions and one electron: two ions are fixed at a distance of 
L = 19.0 ao, the third ion and the electron are free to move in one dimension along the 
line joining the two fixed ions. The moving ion interacts with the fixed ions via a bare 
Coulomb potentials, while the electron-ion interactions are treated as soft-Coulomb po¬ 
tentials, with parameters Rm = 5.0 ao (moving ion), Ri = 3.1 ao (left ion) and Rr = 4.0 oq 
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(right ion). The full Hamiltonian for this system is 


H{r, R) 


1 1 1 1 (^) 
2dr^ 2MdR^^ \^-R\^ \^ + R\ \R - r\ 



(SM.20) 


The nuclear mass is M = 1836. With this choice of parameters the two lowest BO states 
are strongly coupled and there is a weak coupling to the remaining states, as shown in 
Fig. SM.2. The initial wave function is the product of a real-valued normalized Gaussian 
wave packet, centered at Rc = —4.0 oq with variance a = 1/a/2.85 oq, and the second BO 
electronic state. The time-step used for integrating the TDSE is 2.4 x 10“^ fs (or 0.1 a.u.). 
CT-MQC results are obtained by propagating Wraj = 200 trajectories whose initial posi¬ 
tions and momenta are sampled from the Wigner distribution corresponding to the initial 
quantum mechanical density. In this case, the time-step is 1.2 x 10“^ fs (or 0.5 a.u.). The 
forth-order Runge-Kutta algorithm is used to propagate the electronic equation and the 
velocity-Verlet algorithm is used for the nuclear equation. 



R (a.u.) 


Supplemental Material, Figure SM.2. Ground (red) and excited (blue) adiabatic potential energy 
surfaces for the model case (1). The non-adiabatic coupling between the two states is also shown 
(green). 


Test case (2) 

In the Letter we have reported results for the so-called Tully's model #3 [8]. It is a 
two-state model representing an extended coupling region with reflection. This model 
has been well-studied in the literature as it represents a critical test for any new approach 
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to non-adiabatic dynamics based on trajectories. In fact, when a nuclear wave packet 
starts in far negative region and moves with positive but low initial momentum, it is 
partially reflected by the upper surface and partially transmitted on the lower surface. 
This event takes place when the nuclear wave packet reaches the region of space where 
the two adiabatic surfaces start to significantly deviate from each other. The adiabatic 
potential energy surfaces are shown in Fig. SM.3, along with the non-adiabatic coupling. 
The analytic form [8] of the electronic Hamiltonian, in diabatic base, is 


Ffii(i?) = a, H22{R) = 

Hi 2 {R) = 6exp (cR), R < 0, 

HuiR) = b [2 — exp (—ci?)] ,R>0, 
H2i{R) = Hi2{R) 


(SM.21) 


with a = 6 X 10 “"^, 6 = 0 . 1 , c = 0 . 9 . In the Letter we show results for a low initial momen¬ 
tum hko = 10 a.u., which represents a critical situation due to the high probability of the 
reflection chaimel. The initial nuclear wave packet is as Gaussian with width a = 20/ko 
with Rq = — 15.0 a.u. the center of the Gaussian. Exact results are obtained using a wave 
packet propagation scheme where we chose the value 0.1 a.u. for the time-step used to 
integrate the TDSE. GT-MQG results are obtained by propagating a set of Ntraj = 200 tra¬ 
jectories whose initial positions and momenta are sampled from the Wigner distribution 
corresponding to the initial quantum mechanical density. The time-step is 0.5 a.u. for the 
GT-MQG algorithm. 



R (a.u.) 

Supplemental Material, Figure SM.3. Same as in Fig. SM.2 but the model case (2). 


Analysis in momentum space 

We have observed in the Letter that decoherence starts manifesting itself when the nu¬ 
clear density splits in configuration space. The two events, in fact, are strongly related: 
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decoherence cannot be appropriately reproduced if the wave packets propagates "coher¬ 
ently", without splitting. However, the deficiency of TSH in reproducing decoherence 
proves that we can observe splitting without decoherence. Despite this last observation 
we find interesting to look at the splitting of the nuclear density not only in configuration 
space, but also in momentum space. The position representation seems indeed more nat¬ 
ural, due to the fact the BO states are only defined "at fixed nuclear positions" and not 
"at fixed nuclear momenta". In Fig. SM.4 we show the splitting of the nuclear density 
in momentum space, which is perfectly captured by the momentum distribution of the 
classical trajectories. 



-20 0 20 -20 0 20 -20 0 20 
P(a.u.) 


Supplemental Material, Figure SM.4. Splitting of the nuclear density in momentum representation 
(line), compared with quantum-classical results (histogram). 


Computational cost of the CT-MQC scheme 

The actual bottleneck of trajectory-based approaches to non-adiabatic dynamics is the 
computation of electronic structure properties. However, the solution of the full quanfum 
problem requires to treat the nuclei quantum-mechanically and indeed here a trajectory- 
based method is able to cut the computational cost. 

For instance, we have compared in the proton-coupled-electron-transfer problem the 
computational cost needed to solve the TDSE with the CT-MQC equations. The TDSE 
is integrated using split-operator-technique, referenced in the Supplemental Material, 
while the electronic and nuclear equations of the CT-MQC scheme are solved with the 
fourth-order Runge-Kutta algorithm and velocity-Verlet algorithm, respectively. In these 
particular cases, quantum mechanical calculations require a smaller integration time-step 
than the quantum-classical method, i.e. 5 times smaller for the case shown in the Letter. 
Solving the TDSE up to 1000 a.u. with 1000 x 1000 grid points for one nucleus and one 
electron takes approximately 15 minutes on a normal desktop (Intel Core 17-4790 3.60 
GHz), while it takes 1 minute for the CT-MQC with 200 trajectories. To perform the same 


9 
























number of integration steps, the computational cost of the CT-MQC is approximately a 
factor 3 smaller than solving the TDSE. 

As for the scaling of the computational cost in going to higher dimensions, the nu¬ 
merical effort needed to include a larger number of trajectories, in order to appropriately 
sample the configurations space, scales linearly with the number of trajectories. 


TIME-DEPENDENT VECTOR POTENTIAL 

The time-dependent vector potential of the theory is a gauge-dependent quantity Ex¬ 
act quantum-mechanical results are presented in a gauge where the vector potential is 
set to zero (that is possible since the systems under considerations are solved in one- 
dimension). In order to propose a CT-MQC procedure as general as possible, we have 
chosen to work in a different gauge when developing the algorithm, namely 

4'’(*) + 4d(‘) + 53^4' ■ = 0 (SM.22) 

Therefore, as an example, we show in Eig. SM.5 the vector potential calculated within 
the CT-MQC scheme for model case (1). It cannot be compared with exact results, as it 
is identically zero in the quantum case. This is also the reason why the only potential 
shown in the Letter is the gauge-invariant part of the TDPES. The vector potential is non- 


t = 6.1fs t=17.2fs t=24.2fs 



R (ao) R (ao) R (ao) 

Supplemental Material, Figure SM.5. Snapshots of the vector potential taken along the quantum- 
classical evolution. 


zero only after the nuclear density splits in two branches. Notice that it appears as a 
nuclear momentum correction, due to the non-adiabatic nature of the electrons, in the 
nuclear Hamiltonian in Eq. (2) of the Letter. At time t = 24.2 fs in Eig. SM.5, the vector 
potential contributes a negative term only to those trajectories that move on the upper BO 
surface, as its negative part corresponds to the region occupied by the BO-projected wave 
packet "on" the upper BO surface. The opposite situation applies to the positive part of 
the vector potential. This observation is in agreement with previous analysis reported in 
Ref. [3]. 
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